load packages

library(tidyverse)
library(knitr)
devtools::install_github("dcosme/specr", ref = "plotmods")
library(specr)

define palettes

pal_outcome = wesanderson::wes_palette("Zissou1", 3, "continuous")
pal_sca = c(pal_outcome[3], pal_outcome[3], "grey")
pal_condition = wesanderson::wes_palette("Zissou1", 16, "continuous")
pal_condition = c(pal_condition[1], pal_condition[3],
                  pal_condition[14], pal_condition[5],
                  pal_condition[15], pal_condition[16],
                  pal_condition[9], pal_condition[12])

define SCA wrapper functions

plot functions

plot_sca = function(data, combined = TRUE, labels = c("A", "B"), title = FALSE, limits = NULL,
                    point_size = 1, point_alpha = 1, ci_alpha = .5, ci_size = .5, 
                    text_size = 14, title_size = 6, n_breaks = 5,
                    choices  = c("x", "y", "test", "stimulus", "contrast", "controls"),
                    remove_y = FALSE, remove_facet = FALSE) {
  
  if (combined == TRUE) {
      p1 = plot_curve(data, point_size = point_size, point_alpha = point_alpha,
                      ci_alpha = ci_alpha, ci_size = ci_size,
                      limits = limits) +
        geom_hline(aes(yintercept = median(filter(data, y == "body fat percentage")$estimate), linetype = "body fat percentage"), color = pal_outcome[1], size = 1) +
        geom_hline(aes(yintercept = median(filter(data, y == "BMI")$estimate), linetype = "BMI"), color = pal_outcome[2], size = 1) +
        geom_hline(aes(yintercept = median(filter(data, y == "enactment")$estimate), linetype = "enactment"), color = pal_outcome[3], size = 1) +
        scale_linetype_manual(name = "", values = c(2, 2, 2), guide = guide_legend(override.aes = list(color = c(pal_outcome[1], pal_outcome[2], pal_outcome[3])))) +
        labs(x = "", y = "standarized\nregression coefficient\n")  +
        theme(legend.position = "top",
              text = element_text(size = text_size))
      
      if (title == TRUE) {
        if (is.null(limits)) {
          p1 = p1 + annotate("text", -Inf, Inf, label = unique(data$x), fontface = 2, size = title_size,
                       x = 0.5*(1 + nrow(data)), 
                       y = max(data$conf.high))
        } else {
          p1 = p1 + annotate("text", -Inf, Inf, label = unique(data$x), fontface = 2, size = title_size,
                       x = 0.5*(1 + nrow(data)), 
                       y = limits[2])
        }
      }
      
      p2 = plot_choices(data, choices = choices,
                        ignore_vars = c("ROI", "neural signature"), rename_controls = "covariates") +
        labs(x = "\nspecifications (ranked)")  +
        theme(strip.text.x = element_blank(),
              text = element_text(size = text_size))
        
  }
  else {
      p1 = plot_curve(data, point_size = point_size, point_alpha = point_alpha,
                      ci_alpha = ci_alpha, ci_size = ci_size) +
        geom_hline(yintercept = 0, linetype = "solid", color = "black", size = .5) +
        scale_y_continuous(breaks = scales::pretty_breaks(n = n_breaks)) +
        labs(x = "", y = "standarized\nregression coefficient\n") +
        theme(text = element_text(size = text_size))
      
      if (title == TRUE) {
        if (is.null(limits)) {
          p1 = p1 + annotate("text", -Inf, Inf, label = unique(data$y), fontface = 2, size = title_size,
                       x = 0.5*(1 + nrow(data)), 
                       y = max(data$conf.high))
        } else {
          p1 = p1 + annotate("text", -Inf, Inf, label = unique(data$y), fontface = 2, size = title_size,
                       x = 0.5*(1 + nrow(data)), 
                       y = limits[2])
        }
      }
      
      p2 = plot_choices(data, choices = choices,
                        ignore_vars = c("ROI", "neural signature", "not included"), rename_controls = "covariates") +
        labs(x = "\nspecification number (ranked)") +
        theme(strip.text.x = element_blank(),
              text = element_text(size = text_size))
  }
  
  if (remove_y == TRUE) {
    p1 = p1 + labs(y = "")
    
    p2 = p2 + theme(axis.text.y = element_blank(),
                    axis.ticks.y = element_blank()) +
      labs(y = "")
  }

  if (remove_facet == TRUE) {
    p2 = p2 + theme(strip.text.y = element_blank())
    }
  
  plot_specs(plot_a = p1,
             plot_b = p2,
             labels = labels,
             rel_height = c(1, 2))
}

plot_sca_compare = function(data, pointrange = TRUE, labels = c("A", "B"), 
                            rel_heights = c(.75, .25), rel_widths = c(.75, .25), b_limits = NULL,
                            title = FALSE, text_size = 14, title_size = 6, n_rows = 1, n_breaks = 3,
                            remove_x = FALSE, remove_y = FALSE, sig = NULL) {
  
  # define median bootstrapping function
  median_cl_boot = function(x, conf = 0.95, df = TRUE, ci = "low") {
  
    lconf = (1 - conf)/2
    uconf = 1 - lconf
    require(boot)
    bmedian = function(x, ind) median(x[ind])
    bt = boot(x, bmedian, 1000)
    bb = boot.ci(bt, type = "perc")
    
    if (df == TRUE){
      data.frame(y = median(x),
                 ymin = quantile(bt$t, lconf), 
                 ymax = quantile(bt$t, uconf))
      
    } else {
      if (ci == "low") {
        quantile(bt$t, lconf)
      } else {
        quantile(bt$t, uconf)
      }
    }
  }
    
  # merge and tidy for plotting
  plot.data = data %>%
    group_by(x) %>%
    arrange(estimate) %>%
    mutate(specification = row_number()) %>%
    ungroup() %>%
    unique()
  
  # labels
  labs = plot.data %>%
    group_by(x) %>%
    summarize(med = median(estimate),
              low = median_cl_boot(estimate, df = FALSE, ci = "low"),
              high = median_cl_boot(estimate, df = FALSE, ci = "high")) %>%
    mutate(range = max(high) - min(low),
           #estimate = ifelse(med > 0, high + (range / 10), low - (range / 10)),
           estimate = ifelse(med > 0, high + .03, low - .03),
           label = ifelse(x %in% sig, "*", ""))
  
  # plot curves
  if (pointrange == TRUE) {
    a = plot.data %>%
    ggplot(aes(specification, estimate, color = x)) +
      geom_linerange(aes(ymin = conf.low, ymax = conf.high), size = .1) +
      geom_point() +
      geom_hline(yintercept = 0, linetype = "solid", color = "black", size = 1) +
      scale_color_manual(name = "", values = pal_condition) +
      scale_y_continuous(breaks = scales::pretty_breaks(n = n_breaks)) +
      labs(x = "\nspecification number (ranked)", y = "standarized\negression coefficient\n") + 
      theme_minimal() + 
      theme(strip.text = element_blank(), 
            axis.line = element_line("black", size = 0.5), 
            legend.position = c(.5, .1), 
            legend.direction = "horizontal",
            panel.spacing = unit(0.75, "lines"), 
            axis.text = element_text(colour = "black"),
            text = element_text(size = text_size))
    if (title == TRUE) {
      a = a + annotate("text", -Inf, Inf, label = unique(plot.data$y), fontface = 2, size = title_size,
                       x = 0.5*(min(plot.data$specification) + max(plot.data$specification)), 
                       y = max(plot.data$conf.high))
    }
    
  } else {
    a = plot.data %>%
      ggplot(aes(specification, estimate, color = x)) +
      geom_point() +
      geom_hline(yintercept = 0, linetype = "solid", color = "black", size = 1) +
      scale_color_manual(name = "", values = pal_condition) +
      scale_y_continuous(breaks = scales::pretty_breaks(n = n_breaks)) +
      labs(x = "\nspecification number (ranked)", y = "standarized\nregression coefficient\n") + 
      theme_minimal() + 
      theme(strip.text = element_blank(), 
            axis.line = element_line("black", size = 0.5), 
            legend.position = "none", 
            legend.direction = "horizontal",
            panel.spacing = unit(0.75, "lines"), 
            axis.text = element_text(colour = "black"),
            text = element_text(size = text_size))
    if (title == TRUE) {
      a = a + annotate("text", -Inf, Inf, label = unique(plot.data$y), fontface = 2, size = title_size,
                       x = 0.5*(min(plot.data$specification) + max(plot.data$specification)), 
                       y = max(plot.data$estimate))    
      }
  }
  
    b = plot.data %>%
      group_by(x) %>%
      mutate(order = median(estimate)) %>%
      ggplot(aes(reorder(x, order), estimate, fill = x)) +
      stat_summary(fun = "median", geom = "bar") +
      geom_errorbar(data = labs, aes(x = x, ymin = low, ymax = high), width = 0) +
      geom_text(data = labs, aes(label = label, x = x, y = estimate), size = 6) +
      scale_fill_manual(name = "", values = pal_condition) +
      labs(x = "\n", y = "median effect\n") + 
      theme_minimal() + 
      theme(strip.text = element_blank(), 
            axis.line = element_line("black", size = 0.5), 
            legend.position = "none", 
            panel.spacing = unit(0.75, "lines"), 
            axis.text = element_text(colour = "black"),
            text = element_text(size = text_size),
            axis.text.x = element_text(angle = 45, hjust = 1))
    
  if (!is.null(b_limits)) {
    b = b +
      scale_y_continuous(breaks = scales::pretty_breaks(n = n_breaks)) +
      coord_cartesian(ylim = b_limits)
  } else {
    b = b +
      scale_y_continuous(breaks = scales::pretty_breaks(n = n_breaks))
  }
    
  if (n_rows == 1) {
    a = a + theme(legend.position = c(.5, .1))
    b = b + coord_flip() +
      labs(x = "\n", y = "\nmedian effect") + 
      theme(axis.text.x = element_text(angle = 0, hjust = 1),
            axis.text.y = element_blank())
  }     
    

  if (remove_x == TRUE) {
    a = a + labs(x = "")
    
    if (n_rows == 1) {
      b = b + labs(y = "")
    } else {
      b = b + labs(x = "")
    }
  }    
  
  if (remove_y == TRUE) {
    a = a + labs(y = "")
    
    if (n_rows == 1) {
      b = b + labs(x = "")
    } else {
      b = b + labs(y = "")
    }
  }  
    
  cowplot::plot_grid(a, b, labels = labels, rel_heights = rel_heights, rel_widths = rel_widths, nrow = n_rows)
}

summary and bootstapping functions

Process overview

  • use case resampling approach to generate a null dataset for each specification
  • run all model specifications
  • extract the dataset for each model specification
  • force the null on each specification by subtracting the effect of the independent variable of interest (b estimate * x) from the dependent variable (y_value) for each observation in the dataset
  • for each bootstrap iteration, sample with replacement from the null dataset for each specification and run the model specifications to generate a curve
  • extract median estimate, % positive (& significant p < .05), % negative (& significant p < .05)
  • repeat process 1000 times
run_boot = function(data, var, n_shuffles, rerun) {
  
  if (rerun == FALSE & file.exists(sprintf("boot/boot_%s.RDS", var))) {
    out = readRDS(sprintf("boot/boot_%s.RDS", var))
  } else {
    out = sca_boot(data = data, var = var, n_shuffles = n_shuffles)
    saveRDS(out, sprintf("boot/boot_%s.RDS", var))
  }
  return(out)
  
}

sca_boot = function(data, var, n_shuffles, conf.level = 0.95) {
  
  # set seed
  set.seed(63)
  
  # initialize data frame
  df = data.frame()

  # separate uniformity and association data
  data1 = data %>%
    filter(grepl("uniformity|ROI", test)) %>%
    select(-test) %>%
    spread(process, value_std)
  
  data2 = data %>%
    filter(grepl("association|ROI", test)) %>%
    select(-test) %>%
    spread(process, value_std)
  
  # define control variables, subsets, and run scas
  # define controls
  if (length(var) > 1) {
    controls = c(unique(data$process), "sex", "restraint")
  } else {
    controls = c(unique(data$process)[!unique(data$process) %in% var], "sex", "restraint")
  }
  
  # define subsets
  subsets = list(stimulus = unique(data$stimulus),
                 contrast = unique(data$contrast))
  
  # run enact_prop models separately
  if (any(outcome == "enact_prop" & length(outcome) > 1)) {
    # define controls
    controls_enact = controls[controls != "sex"]
    
    # run scas
    results_data1_enact = run_specs(df = data1,
                              y = "enact_prop", 
                              x = var,
                              controls = controls_enact,
                              model = c("lm"), 
                              subsets = subsets,
                              keep.results = TRUE,
                              keep.formula = TRUE) %>%
      mutate(test = ifelse(grepl("ROI", controls), "ROI", 
                    ifelse(grepl("ROI", x) & grepl("no covariates", controls), "ROI", 
                    ifelse(controls == "craving_regulation_PEV", "neural signature",     
                    ifelse(grepl("craving_regulation", x) & grepl("no covariates", controls), "neural signature", "uniformity")))))
    
    results_data2_enact = run_specs(df = data2,
                              y = "enact_prop", 
                              x = var,
                              controls = controls_enact,
                              model = c("lm"), 
                              subsets = subsets,
                              keep.results = TRUE,
                              keep.formula = TRUE)  %>%
      mutate(test = ifelse(grepl("ROI", controls), "ROI", 
                    ifelse(grepl("ROI", x) & grepl("no covariates", controls), "ROI", 
                    ifelse(controls == "craving_regulation_PEV", "neural signature",     
                    ifelse(grepl("craving_regulation", x) & grepl("no covariates", controls), "neural signature", "association")))))
    
  }
  
    if (any(outcome == "enact_prop") & length(outcome) == 1) {
      
      # define controls
      controls = controls[controls != "sex"]
     
      # run scas
      results_data1 = run_specs(df = data1,
                                y = outcome, 
                                x = var,
                                controls = controls,
                                model = c("lm"), 
                                subsets = subsets,
                                keep.results = TRUE,
                                keep.formula = TRUE)  %>%
       mutate(test = ifelse(grepl("ROI", controls), "ROI", 
                     ifelse(grepl("ROI", x) & grepl("no covariates", controls), "ROI", 
                     ifelse(controls == "craving_regulation_PEV", "neural signature",     
                     ifelse(grepl("craving_regulation", x) & grepl("no covariates", controls), "neural signature", "uniformity")))))
      
      results_data2 = run_specs(df = data2,
                                y = outcome, 
                                x = var,
                                controls = controls,
                                model = c("lm"), 
                                subsets = subsets,
                                keep.results = TRUE,
                                keep.formula = TRUE)  %>%
      mutate(test = ifelse(grepl("ROI", controls), "ROI", 
                    ifelse(grepl("ROI", x) & grepl("no covariates", controls), "ROI", 
                    ifelse(controls == "craving_regulation_PEV", "neural signature",     
                    ifelse(grepl("craving_regulation", x) & grepl("no covariates", controls), "neural signature", "association")))))
    
    } else {
      
      # run scas
      results_data1 = run_specs(df = data1,
                             y = outcome[!grepl("enact_prop", outcome)], 
                             x = var,
                             controls = controls,
                             model = c("lm"), 
                             subsets = subsets,
                             keep.results = TRUE,
                             keep.formula = TRUE) %>%
      mutate(test = ifelse(grepl("ROI", controls), "ROI", 
                   ifelse(grepl("ROI", x) & grepl("no covariates", controls), "ROI", 
                   ifelse(controls == "craving_regulation_PEV", "neural signature",     
                   ifelse(grepl("craving_regulation", x) & grepl("no covariates", controls), "neural signature", "uniformity")))))
      
      results_data2 = run_specs(df = data2,
                             y = outcome[!grepl("enact_prop", outcome)], 
                             x = var,
                             controls = controls,
                             model = c("lm"), 
                             subsets = subsets,
                             keep.results = TRUE,
                             keep.formula = TRUE)  %>%
      mutate(test = ifelse(grepl("ROI", controls), "ROI", 
                   ifelse(grepl("ROI", x) & grepl("no covariates", controls), "ROI", 
                   ifelse(controls == "craving_regulation_PEV", "neural signature",     
                   ifelse(grepl("craving_regulation", x) & grepl("no covariates", controls), "neural signature", "association")))))
    }
    
  # merge results dataframes
  merged = bind_rows(results_data1, results_data2)
  
  if (any(outcome == "enact_prop" & length(outcome) > 1)) {
    merged = bind_rows(merged, results_data1_enact, results_data2_enact)
  }
    
  # create null dataset by subtracting the effect of x (estimate * x) from the dependent variable (y_value)
  null_data = merged %>%
    filter(!x == controls) %>%
    mutate(stimulus = ifelse(grepl("snack", subsets), "snack",
                      ifelse(grepl("meal", subsets), "meal",
                      ifelse(grepl("dessert", subsets), "dessert",
                      ifelse(grepl("food", subsets), "food (average)", "all")))),
           contrast = ifelse(grepl("nature", subsets), "nature",
                      ifelse(grepl("rest", subsets), "rest",
                      ifelse(grepl("average", subsets), "average", "all"))),
           duplicated = ifelse(duplicated(estimate) & duplicated(std.error), 1, 0)) %>%
    filter(!duplicated) %>%
    unique() %>%
    ungroup() %>%
    filter(!stimulus == "all" & !contrast == "all") %>%
    rownames_to_column() %>%
    rename("model_number" = rowname) %>%
    group_by(model_number) %>%
    mutate(data = list(res[[1]][["model"]])) %>%
    select(-res) %>%
    unnest() %>%
    group_by(model_number) %>%
    mutate(obs_number = row_number()) %>%
    ungroup() %>%
    gather(outcome, y_value, bmi, fat, enact_prop) %>%
    mutate(y_null = y_value - (estimate * !!as.name(var))) %>%
    select(-y_value, -estimate) %>%
    spread(outcome, y_null) %>%
    select(-c(std.error, statistic, p.value, conf.low, conf.high, obs, obs_number, test))
  
  for (i in 1:n_shuffles){
    
    sampled = null_data %>%
      group_by(model_number) %>%
      mutate(bmi = ifelse(y == "bmi", sample(bmi, replace = TRUE), NA),
             fat = ifelse(y == "fat", sample(fat, replace = TRUE), NA),
             enact_prop = ifelse(y == "enact_prop", sample(enact_prop, replace = TRUE), NA)) %>%
      group_by(model_number, x, y, model, formula, controls, random_effects, subsets, stimulus, contrast) %>%
      nest()
    
    results = sampled %>%
      mutate(res = map2(.x = .data$data,
                        .y = .data$formula,
                        .f = ~ lm(.y, data = .x)),
             coefs = map(.data$res,
                         broom::tidy,
                         conf.int = TRUE,
                         conf.level = conf.level),
             obs = map(.data$res, nobs)) %>%
      unnest(.data$coefs) %>%
      unnest(.data$obs) %>%
      filter(.data$term == .data$x) %>%
      select(-.data$term, -.data$res, -.data$formula) %>%
      mutate(pos = ifelse(estimate > 0, 1, 0),
             neg = ifelse(estimate < 0, 1, 0),
             sig = ifelse(p.value < .05, 1, 0),
             pos_sig = ifelse(pos == 1 & sig == 1, 1, 0),
             neg_sig = ifelse(neg == 1 & sig == 1, 1, 0)) %>%
      group_by(x, y) %>%
      summarize(median = median(estimate),
                n = n(),
                n_positive = sum(pos),
                n_negative = sum(neg),
                n_significant = sum(sig),
                n_positive_sig = sum(pos_sig),
                n_negative_sig = sum(neg_sig)) %>%
      mutate(shuffle = i) %>%
      select(shuffle, everything()) %>%
      ungroup() %>%
      mutate_if(is.character, stringr::str_replace_all, pattern = "enact_prop", replacement = "enactment") %>%
      mutate_if(is.character, stringr::str_replace_all, pattern = "bmi", replacement = "BMI") %>%
      mutate_if(is.character, stringr::str_replace_all, pattern = "fat", replacement = "body fat percentage") %>%
      mutate_if(is.character, stringr::str_replace_all, pattern = "_", replacement = " ")
  
    # join to data frame
    df = rbind(df, as.data.frame(results))
    rm(results)
    
  }
  
  return(df)
  
}
      
run_sca = function(data, var, outcome, boot = FALSE, n_shuffles = 1, rerun = FALSE) {
    
    # set seed
    set.seed(63)
  
    # define median bootstrapping function
    median_cl_boot = function(estimate, conf = 0.95) {
  
      lconf = (1 - conf)/2
      uconf = 1 - lconf
      require(boot)
      bmedian = function(estimate, ind) median(estimate[ind])
      bt = boot(estimate, bmedian, 1000)
      bb = boot.ci(bt, type = "perc")
      data.frame(obs_median = median(estimate), 
                 obs_conf_low = quantile(bt$t, lconf), 
                 obs_conf_high = quantile(bt$t, uconf))
    
    }
    
    # separate uniformity and association data
    data1 = data %>%
      filter(grepl("uniformity|ROI", test)) %>%
      select(-test) %>%
      spread(process, value_std)
    
    data2 = data %>%
      filter(grepl("association|ROI", test)) %>%
      select(-test) %>%
      spread(process, value_std)
    
    # define control variables, subsets, and run scas
    # define controls
    if (length(var) > 1) {
      controls = c(unique(data$process), "sex", "restraint")
    } else {
      controls = c(unique(data$process)[!unique(data$process) %in% var], "sex", "restraint")
    }
    
    # define subsets
    subsets = list(stimulus = unique(data$stimulus),
                   contrast = unique(data$contrast))
  
    # run enact_prop models separately
    if (any(outcome == "enact_prop") & length(outcome) > 1) {
      
      # define controls
      controls_enact = controls[controls != "sex"]
      
      # run scas
      results_data1_enact = run_specs(df = data1,
                                y = "enact_prop", 
                                x = var,
                                controls = controls_enact,
                                model = c("lm"), 
                                subsets = subsets) %>%
        mutate(test = ifelse(grepl("ROI", controls), "ROI", 
                      ifelse(grepl("ROI", x) & grepl("no covariates", controls), "ROI", 
                      ifelse(controls == "craving_regulation_PEV", "neural signature",     
                      ifelse(grepl("craving_regulation", x) & grepl("no covariates", controls), "neural signature", "uniformity")))))
      
      results_data2_enact = run_specs(df = data2,
                                y = "enact_prop", 
                                x = var,
                                controls = controls_enact,
                                model = c("lm"), 
                                subsets = subsets)  %>%
        mutate(test = ifelse(grepl("ROI", controls), "ROI", 
                      ifelse(grepl("ROI", x) & grepl("no covariates", controls), "ROI", 
                      ifelse(controls == "craving_regulation_PEV", "neural signature",     
                      ifelse(grepl("craving_regulation", x) & grepl("no covariates", controls), "neural signature", "association")))))
      
    }
    
    if (any(outcome == "enact_prop") & length(outcome) == 1) {
      
      # define controls
      controls = controls[controls != "sex"]
     
      # run scas
      results_data1 = run_specs(df = data1,
                                y = outcome, 
                                x = var,
                                controls = controls,
                                model = c("lm"), 
                                subsets = subsets)  %>%
       mutate(test = ifelse(grepl("ROI", controls), "ROI", 
                     ifelse(grepl("ROI", x) & grepl("no covariates", controls), "ROI", 
                     ifelse(controls == "craving_regulation_PEV", "neural signature",     
                     ifelse(grepl("craving_regulation", x) & grepl("no covariates", controls), "neural signature", "uniformity")))))
      
      results_data2 = run_specs(df = data2,
                                y = outcome, 
                                x = var,
                                controls = controls,
                                model = c("lm"), 
                                subsets = subsets)  %>%
      mutate(test = ifelse(grepl("ROI", controls), "ROI", 
                    ifelse(grepl("ROI", x) & grepl("no covariates", controls), "ROI", 
                    ifelse(controls == "craving_regulation_PEV", "neural signature",     
                    ifelse(grepl("craving_regulation", x) & grepl("no covariates", controls), "neural signature", "association")))))
    
    } else {
      
      # run scas
      results_data1 = run_specs(df = data1,
                             y = outcome[!grepl("enact_prop", outcome)], 
                             x = var,
                             controls = controls,
                             model = c("lm"), 
                             subsets = subsets) %>%
      mutate(test = ifelse(grepl("ROI", controls), "ROI", 
                   ifelse(grepl("ROI", x) & grepl("no covariates", controls), "ROI", 
                   ifelse(controls == "craving_regulation_PEV", "neural signature",     
                   ifelse(grepl("craving_regulation", x) & grepl("no covariates", controls), "neural signature", "uniformity")))))
      
      results_data2 = run_specs(df = data2,
                             y = outcome[!grepl("enact_prop", outcome)], 
                             x = var,
                             controls = controls,
                             model = c("lm"), 
                             subsets = subsets)  %>%
      mutate(test = ifelse(grepl("ROI", controls), "ROI", 
                   ifelse(grepl("ROI", x) & grepl("no covariates", controls), "ROI", 
                   ifelse(controls == "craving_regulation_PEV", "neural signature",     
                   ifelse(grepl("craving_regulation", x) & grepl("no covariates", controls), "neural signature", "association")))))
    }
    
    # merge results dataframes
    merged = bind_rows(results_data1, results_data2)
    
    if (any(outcome == "enact_prop" & length(outcome) > 1)) {
      merged = bind_rows(merged, results_data1_enact, results_data2_enact)
    }
    
    # extract info & exclude duplicates
    results = merged %>%
      filter(!x == controls) %>%
      mutate(stimulus = ifelse(grepl("snack", subsets), "snack",
                        ifelse(grepl("meal", subsets), "meal",
                        ifelse(grepl("dessert", subsets), "dessert",
                        ifelse(grepl("food", subsets), "food (average)", "all")))),
             contrast = ifelse(grepl("nature", subsets), "nature",
                        ifelse(grepl("rest", subsets), "rest",
                        ifelse(grepl("average", subsets), "average", "all"))),
             duplicated = ifelse(duplicated(estimate) & duplicated(std.error), 1, 0)) %>%
      filter(!duplicated) %>%
      unique() %>%
      ungroup() %>%
      filter(!stimulus == "all" & !contrast == "all") %>%
      mutate_if(is.character, stringr::str_replace_all, pattern = "enact_prop", replacement = "enactment") %>%
      mutate_if(is.character, stringr::str_replace_all, pattern = "bmi", replacement = "BMI") %>%
      mutate_if(is.character, stringr::str_replace_all, pattern = "fat", replacement = "body fat percentage") %>%
      mutate_if(is.character, stringr::str_replace_all, pattern = "_", replacement = " ")
    
    median_ci = results %>%
      group_by(x, y) %>%
      do({
        median_cl_boot(.$estimate)
      })
        
    summary = results %>%
      mutate(pos = ifelse(estimate > 0, 1, 0),
             neg = ifelse(estimate < 0, 1, 0),
             sig = ifelse(p.value < .05, 1, 0),
             pos_sig = ifelse(pos == 1 & sig == 1, 1, 0),
             neg_sig = ifelse(neg == 1 & sig == 1, 1, 0)) %>%
      group_by(x, y) %>%
      summarize(obs_n = n(),
                obs_n_positive = sum(pos),
                obs_n_negative = sum(neg),
                obs_n_significant = sum(sig),
                obs_n_positive_sig = sum(pos_sig),
                obs_n_negative_sig = sum(neg_sig)) %>%
      left_join(., median_ci, by = c("x", "y")) %>%
      select(x, y, obs_median, obs_conf_low, obs_conf_high, everything())
      

  
    if (boot == TRUE) {
      boot_out = run_boot(data, var, n_shuffles, rerun)
      return(list(results = results, summary = summary, boot = boot_out))
    } else {
      return(list(results = results, summary = summary))
    }
    
}

boot_stats = function(data) {
  
  data$boot %>%
    mutate(y = ifelse(grepl("%", y), "body fat percentage", y)) %>%
    left_join(., data$summary, c("x", "y")) %>%
    mutate(extreme_median = ifelse(obs_median > 0 & median >= obs_median, 1,
                            ifelse(obs_median < 0 & median <= obs_median, 1, 0)),
           extreme_positive = ifelse(n_positive >= obs_n_positive, 1, 0),
           extreme_negative = ifelse(n_negative >= obs_n_negative, 1, 0),
           extreme_positive_sig = ifelse(n_positive_sig >= obs_n_positive_sig, 1, 0),
           extreme_negative_sig = ifelse(n_negative_sig >= obs_n_negative_sig, 1, 0)) %>%
    group_by(x, y) %>%
    summarize(n = n(),
              extreme_median = sum(extreme_median),
              extreme_positive = sum(extreme_positive),
              extreme_negative = sum(extreme_negative),
              extreme_positive_sig = sum(extreme_positive_sig),
              extreme_negative_sig = sum(extreme_negative_sig)) %>%
    gather(variable, n_extreme, contains("extreme")) %>%
    mutate(p_value = n_extreme / n,
           p_value = ifelse(p_value == 1.000, "1.000", 
               ifelse(p_value < .001, "< .001", gsub("0.(.*)", ".\\1", sprintf("%.3f", p_value))))) %>%
    select(x, y, variable, p_value) %>%
    spread(variable, p_value) %>%
    left_join(., select(ungroup(data$summary), x, y, obs_median, obs_conf_low, obs_conf_high, obs_n_positive_sig, obs_n_negative_sig), by = c("x", "y")) %>%
    mutate(`Mdn [95% CI]` = sprintf("%.2f [%.2f, %.2f]", obs_median, obs_conf_low, obs_conf_high)) %>%
    select(y, `Mdn [95% CI]`, extreme_median, obs_n_positive_sig, extreme_positive_sig, obs_n_negative_sig, extreme_negative_sig) %>%
    rename("Mdn p" = extreme_median,
           "Positive N" = obs_n_positive_sig,
           "Positive p" = extreme_positive_sig,
           "Negative N" = obs_n_negative_sig,
           "Negative p" = extreme_negative_sig) %>%
    gather(var, val, -c(x, y)) %>%
    unite(y, y, var, sep = " ") %>%
    spread(y, val)
}

load and tidy data

Read betas_dataset.RDS and dots_dataset.RDS saved from the xval_models.Rmd

  • Ran uniformity and association test data separately because the models were breaking when trying to include test type as a covariate
  • Code sex as -.5 (female), .5 (male)
source("load_data.R")

ind_diffs1 = ind_diffs %>%
  mutate(bmi = as.numeric(scale(bmi)),
         fat = as.numeric(scale(fat)),
         restraint = as.numeric(scale(restraint)),
         sex = ifelse(gender == 1, -.5,
                  ifelse(gender == 2, .5, gender))) %>%
  select(-gender)

ema_enact1 = ema_enact %>%
  mutate(enact_prop = as.numeric(scale(enact_prop)))

betas_std = readRDS("betas_dataset.RDS") %>%
  filter(session == "all") %>%
  group_by(subjectID, process, condition, control) %>%
  mutate(value_std = mean(meanPE_std, na.rm = TRUE)) %>%
  ungroup() %>%
  select(subjectID, process, condition, control, value_std) %>%
  unique() %>%
  mutate(process = sprintf("%s_ROI", process)) %>%
  filter(grepl("snack|meal|dessert|food", condition)) %>%
  spread(control, value_std) %>%
  mutate(average = (nature + rest) / 2) %>%
  gather(control, value_std, nature, rest, average)

data_unif = readRDS("dots_dataset.RDS") %>%
  filter(session == "all" & mask == "unmasked") %>%
  filter((test == "uniformity" & !process == "craving_regulation") | 
           (test == "association" & process == "craving_regulation"))  %>%
  mutate(process = sprintf("%s_PEV", process)) %>%
  rename("value_std" = dotProduct_std) %>%
  select(subjectID, process, condition, control, value_std) %>%
  unique() %>%
  filter(grepl("snack|meal|dessert|food", condition)) %>%
  spread(control, value_std) %>%
  mutate(average = (nature + rest) / 2) %>%
  gather(control, value_std, nature, rest, average) %>%
  bind_rows(betas_std) %>%
  left_join(., ind_diffs1, by = "subjectID") %>%
  left_join(., ema_enact1, by = "subjectID") %>%
  select(-c(sample, DBIC_ID, age)) %>% #gender
  rename("contrast" = control,
         "stimulus" = condition) %>%
  mutate(stimulus = as.factor(stimulus),
         contrast = as.factor(contrast)) %>%
  spread(process, value_std) %>%
  mutate(stimulus = gsub("food", "food (average)", stimulus))

data_assoc = readRDS("dots_dataset.RDS") %>%
  filter(session == "all" & mask == "unmasked" & test == "association") %>%
  mutate(process = sprintf("%s_PEV", process)) %>%
  rename("value_std" = dotProduct_std) %>%
  select(subjectID, process, condition, control, value_std) %>%
  unique() %>%
  filter(grepl("snack|meal|dessert|food", condition)) %>%
  spread(control, value_std) %>%
  mutate(average = (nature + rest) / 2) %>%
  gather(control, value_std, nature, rest, average) %>%
  bind_rows(betas_std) %>%
  left_join(., ind_diffs1, by = "subjectID") %>%
  left_join(., ema_enact1, by = "subjectID") %>%
  select(-c(sample, DBIC_ID, age)) %>% #gender
  rename("contrast" = control,
         "stimulus" = condition) %>%
  mutate(stimulus = as.factor(stimulus),
         contrast = as.factor(contrast)) %>%
  spread(process, value_std) %>%
  mutate(stimulus = gsub("food", "food (average)", stimulus))

data_unif_gathered = data_unif %>%
  gather(process, value_std, contains("ROI"), contains("PEV")) %>%
  mutate(test = ifelse(grepl("PEV", process), "uniformity", "ROI"))

data_assoc_gathered = data_assoc %>%
  gather(process, value_std, contains("ROI"), contains("PEV")) %>%
  mutate(test = ifelse(grepl("PEV", process), "association", "ROI"))

data_combined = bind_rows(data_unif_gathered, data_assoc_gathered) %>%
  unique()

head(data_combined)

large coefficient models

specify modeling approach

  • Also include the following subsets: stimulus type ( meal, snack, dessert, or food (average)) and contrast (nature, rest, or)
setup_specs(y = c("bmi", "fat", "enact_prop"),
            x = c("reward_ROI", "reward_PEV",
                  "value_ROI", "value_PEV",
                  "cognitive_control_ROI", "cognitive_control_PEV",
                  "craving_PEV", "craving_regulation_PEV"),
            controls = c("reward_ROI", "reward_PEV",
                         "value_ROI", "value_PEV",
                         "cognitive_control_ROI", "cognitive_control_PEV",
                         "craving_PEV", "craving_regulation_PEV"),
            model = c("lm")) %>%
  filter(!x == controls)

define common variables

data = data_combined
var = c("reward_ROI", "reward_PEV",
        "value_ROI", "value_PEV",
        "cognitive_control_ROI", "cognitive_control_PEV",
        "craving_PEV", "craving_regulation_PEV")

run SCAs

BMI

# define outcome
outcome = "bmi"

# run sca
output = run_sca(data = data, var = var, outcome = outcome, boot = FALSE)

bmi_output_plot = output$results %>%
  mutate(`a priori` = ifelse(x == "value ROI" & controls == "reward ROI" &
                               stimulus == "food (average)" & contrast == "nature", "ROI model", 
                      ifelse(x == "reward ROI" & controls == "value ROI" &
                               stimulus == "food (average)" & contrast == "nature", "ROI model", 
                      ifelse(x == "craving regulation PEV" & controls == "no covariates" & 
                               stimulus == "food (average)" & contrast == "nature", "PEV model", "not included"))))
  
# plot
plot_sca(data = bmi_output_plot, combined = FALSE, choices = c("x", "test", "stimulus", "contrast", "controls", "a priori"))

(bmi_compare = plot_sca_compare(data = bmi_output_plot, pointrange = FALSE))

plot_decisiontree(bmi_output_plot, legend = TRUE)

% fat

# define outcome
outcome = "fat"

# run sca
output = run_sca(data = data, var = var, outcome = outcome, boot = FALSE, n_shuffles = n_shuffles, rerun = FALSE)

fat_output_plot = output$results %>%
  mutate(`a priori` = ifelse(x == "reward ROI" & controls == "no covariates" &
                               stimulus == "food (average)" & contrast == "nature", "ROI model", 
                      ifelse(x == "craving PEV" & controls == "no covariates" & 
                               stimulus == "food (average)" & contrast == "nature" & test == "association", "PEV model", "not included")),
         y = ifelse(grepl("%", y), "body fat percentage", y))

# plot
plot_sca(data = fat_output_plot, combined = FALSE, choices = c("x", "test", "stimulus", "contrast", "controls", "a priori"))

(fat_compare = plot_sca_compare(data = fat_output_plot, pointrange = FALSE))

plot_decisiontree(fat_output_plot, legend = TRUE)

enactment

# define outcome
outcome = "enact_prop"

# run sca
output = run_sca(data = data, var = var, outcome = outcome, boot = FALSE, n_shuffles = n_shuffles, rerun = FALSE)

enact_output_plot = output$results %>%
  mutate(`a priori` = ifelse(x == "value ROI" & controls == "no covariates" &
                               stimulus == "food (average)" & contrast == "nature", "ROI model", 
                      ifelse(x == "reward PEV" & controls == "no covariates" & 
                               stimulus == "food (average)" & contrast == "nature" & test == "association", "PEV model", "not included")))

# plot
plot_sca(data = enact_output_plot, combined = FALSE, choices = c("x", "test", "stimulus", "contrast", "controls", "a priori"))

(enact_compare = plot_sca_compare(data = enact_output_plot, pointrange = FALSE))

plot_decisiontree(enact_output_plot, legend = TRUE)

combined plots

bmi_sca = plot_sca(data = bmi_output_plot, combined = FALSE, labels = c("A", "B"), title = TRUE,
                   point_size = .2, point_alpha = 1, ci_size = .1, ci_alpha = .5, text_size = 18,
                   remove_facet = TRUE, n_breaks = 5,
                   choices = c("x", "test", "stimulus", "contrast", "a priori"))
fat_sca = plot_sca(data = fat_output_plot, combined = FALSE, labels = c("C", "D"), title = TRUE,
                   point_size = .2, point_alpha = 1, ci_size = .1, ci_alpha = .5, text_size = 18,
                   remove_y = TRUE, remove_facet = TRUE, n_breaks = 5,
                   choices = c("x", "test", "stimulus", "contrast", "a priori"))
enact_sca = plot_sca(data = enact_output_plot, combined = FALSE, labels = c("E", "F"), title = TRUE,
                     point_size = .2, point_alpha = 1, ci_size = .1, ci_alpha = .5, text_size = 18,
                     remove_y = TRUE, n_breaks = 5,
                     choices = c("x", "test", "stimulus", "contrast", "a priori"))

cowplot::plot_grid(bmi_sca, fat_sca, enact_sca, ncol = 3, rel_widths = c(.38, .31, .31))

bmi_compare = plot_sca_compare(data = bmi_output_plot, pointrange = FALSE, title = TRUE,
                               labels = c("A", "B"), rel_heights = c(.5, .5), text_size = 18,
                               n_rows = 2, n_breaks = 3,b_limits   = c(-.25, .3),
                               sig = c("reward ROI",
                                       "value ROI",
                                       "cognitive control PEV",
                                       "craving PEV", "craving regulation PEV"))
fat_compare = plot_sca_compare(data = fat_output_plot, pointrange = FALSE, title = TRUE,
                               labels = c("C", "D"), rel_heights = c(.5, .5), text_size = 18,
                               remove_y = TRUE, n_rows = 2, n_breaks = 3, b_limits  = c(-.25, .3),
                               sig = c("value ROI", 
                                       "cognitive control ROI",
                                       "craving PEV", "craving regulation PEV"))
enact_compare = plot_sca_compare(data = enact_output_plot, pointrange = FALSE, title = TRUE,
                                 labels = c("E", "F"), rel_heights = c(.5, .5), text_size = 18,
                                 remove_y = TRUE, n_rows = 2, n_breaks = 3, b_limits = c(-.25, .3),
                                 sig = c("reward ROI", "reward PEV",
                                       "value ROI", "value PEV",
                                       "cognitive control ROI",
                                       "craving PEV"))

cowplot::plot_grid(bmi_compare, fat_compare, enact_compare, ncol = 3)

neural covariate models

define common variables

data = data_combined
n_shuffles = 1000
outcome = c("bmi", "fat", "enact_prop")

run SCAs

reward ROI

# define variable
var = "reward_ROI"

# run SCA
reward_ROI_output = run_sca(data = data, var = var, outcome = outcome, boot = TRUE, n_shuffles = n_shuffles, rerun = FALSE)

reward_ROI_output$summary 
# stats
(reward_ROI = boot_stats(data = reward_ROI_output))
# plot
plot_sca(data = reward_ROI_output$results, combined = TRUE, title = TRUE)

reward PEV

# define variable
var = "reward_PEV"

# run SCA
reward_PEV_output = run_sca(data = data, var = var, outcome = outcome, boot = TRUE, n_shuffles = n_shuffles, rerun = FALSE)
reward_PEV_output$summary 
# stats
(reward_PEV = boot_stats(data = reward_PEV_output))
# plot
plot_sca(data = reward_PEV_output$results, combined = TRUE, title = TRUE)

value ROI

# define variable
var = "value_ROI"

# run SCA
value_ROI_output = run_sca(data = data, var = var, outcome = outcome, boot = TRUE, n_shuffles = n_shuffles, rerun = FALSE)
value_ROI_output$summary 
# stats
(value_ROI = boot_stats(data = value_ROI_output))
# plot
plot_sca(data = value_ROI_output$results, combined = TRUE, title = TRUE)

value PEV

# define variable
var = "value_PEV"

# run SCA
value_PEV_output = run_sca(data = data, var = var, outcome = outcome, boot = TRUE, n_shuffles = n_shuffles, rerun = FALSE)
value_PEV_output$summary 
# stats
(value_PEV = boot_stats(data = value_PEV_output))
# plot
plot_sca(data = value_PEV_output$results, combined = TRUE, title = TRUE)

cognitive control ROI

# define variable
var = "cognitive_control_ROI"

# run SCA
cognitive_control_ROI_output = run_sca(data = data, var = var, outcome = outcome, boot = TRUE, n_shuffles = n_shuffles, rerun = FALSE)
cognitive_control_ROI_output$summary 
# stats
(cognitive_control_ROI = boot_stats(data = cognitive_control_ROI_output))
# plot
plot_sca(data = cognitive_control_ROI_output$results, combined = TRUE, title = TRUE)

cognitive control PEV

# define variable
var = "cognitive_control_PEV"

# run SCA
cognitive_control_PEV_output = run_sca(data = data, var = var, outcome = outcome, boot = TRUE, n_shuffles = n_shuffles, rerun = FALSE)
cognitive_control_PEV_output$summary 
# stats
(cognitive_control_PEV = boot_stats(data = cognitive_control_PEV_output))
# plot
plot_sca(data = cognitive_control_PEV_output$results, combined = TRUE, title = TRUE)

craving PEV

# define variable
var = "craving_PEV"

# run SCA
craving_PEV_output = run_sca(data = data, var = var, outcome = outcome, boot = TRUE, n_shuffles = n_shuffles, rerun = FALSE)
craving_PEV_output$summary 
# stats
(craving_PEV = boot_stats(data = craving_PEV_output))
# plot
plot_sca(data = craving_PEV_output$results, combined = TRUE, title = TRUE)

craving regulation PEV

# define variable
var = "craving_regulation_PEV"

# run SCA
craving_regulation_PEV_output = run_sca(data = data, var = var, outcome = outcome, boot = TRUE, n_shuffles = n_shuffles, rerun = FALSE)
craving_regulation_PEV_output$summary 
# stats
(craving_regulation_PEV = boot_stats(data = craving_regulation_PEV_output))
# plot
plot_sca(data = craving_regulation_PEV_output$results, combined = TRUE, title = TRUE)

plot combined SCAs

reward

limits = c(min(reward_PEV_output$results$conf.low), max(reward_PEV_output$results$conf.high))

a = plot_sca(data = mutate(reward_ROI_output$results, 
                           controls = ifelse(controls == "reward PEV", "reward PEV / ROI", controls)), 
             combined = TRUE, title = TRUE, labels = c("A", "B"),
             choices = c("y", "test", "stimulus", "contrast", "controls"),
             remove_facet = TRUE, text_size = 18, limits = limits)

b = plot_sca(data = reward_PEV_output$results, combined = TRUE, title = TRUE, labels = c("C", "D"),
         choices = c("y", "test", "stimulus", "contrast", "controls"),
         remove_y = TRUE, text_size = 18, limits = limits)

cowplot::plot_grid(a, b, ncol = 2, rel_widths = c(.54, .46))

value

limits = c(min(value_PEV_output$results$conf.low), max(value_PEV_output$results$conf.high))

a = plot_sca(data = mutate(value_ROI_output$results, 
                           controls = ifelse(controls == "value PEV", "value PEV / ROI", controls)), 
             combined = TRUE, title = TRUE, labels = c("A", "B"),
             choices = c("y", "test", "stimulus", "contrast", "controls"),
             remove_facet = TRUE, text_size = 18, limits = limits)

b = plot_sca(data = value_PEV_output$results, combined = TRUE, title = TRUE, labels = c("C", "D"),
         choices = c("y", "test", "stimulus", "contrast", "controls"),
         remove_y = TRUE, text_size = 18, limits = limits)

cowplot::plot_grid(a, b, ncol = 2, rel_widths = c(.53, .47))

cognitive control

limits = c(min(cognitive_control_ROI_output$results$conf.low), max(cognitive_control_PEV_output$results$conf.high))

a = plot_sca(data = mutate(cognitive_control_ROI_output$results, 
                           controls = ifelse(controls == "cognitive control PEV", "cognitive control PEV / ROI", controls)), 
             combined = TRUE, title = TRUE, labels = c("A", "B"),
             choices = c("y", "test", "stimulus", "contrast", "controls"),
             remove_facet = TRUE, text_size = 18, limits = limits)

b = plot_sca(data = cognitive_control_PEV_output$results, combined = TRUE, title = TRUE, labels = c("C", "D"),
         choices = c("y", "test", "stimulus", "contrast", "controls"),
         remove_y = TRUE, text_size = 18, limits = limits)

cowplot::plot_grid(a, b, ncol = 2, rel_widths = c(.57, .43))

craving & craving regulation

limits = c(min(craving_PEV_output$results$conf.low), max(craving_PEV_output$results$conf.high))

a = plot_sca(data = mutate(craving_PEV_output$results, 
                           controls = ifelse(controls == "craving regulation PEV", "craving regulation / craving PEV", controls)), 
             combined = TRUE, title = TRUE, labels = c("A", "B"),
             choices = c("y", "test", "stimulus", "contrast", "controls"),
             remove_facet = TRUE, text_size = 18, limits = limits)

b = plot_sca(data = craving_regulation_PEV_output$results, combined = TRUE, title = TRUE, labels = c("C", "D"),
         choices = c("y", "test", "stimulus", "contrast", "controls"),
         remove_y = TRUE, text_size = 18, limits = limits)

cowplot::plot_grid(a, b, ncol = 2, rel_widths = c(.57, .43))

combined table

bind_rows(reward_ROI, reward_PEV, craving_PEV, cognitive_control_ROI, cognitive_control_PEV, craving_regulation_PEV, value_ROI, value_PEV) %>%
  kable(format = "pandoc", digits = 2)
x BMI Mdn [95% CI] BMI Mdn p BMI Negative N BMI Negative p BMI Positive N BMI Positive p body fat percentage Mdn [95% CI] body fat percentage Mdn p body fat percentage Negative N body fat percentage Negative p body fat percentage Positive N body fat percentage Positive p enactment Mdn [95% CI] enactment Mdn p enactment Negative N enactment Negative p enactment Positive N enactment Positive p
reward ROI -0.13 [-0.14, -0.11] < .001 22 < .001 0 1.000 -0.02 [-0.03, -0.01] .012 0 1.000 0 1.000 0.19 [0.17, 0.21] < .001 0 1.000 24 < .001
reward PEV -0.00 [-0.01, 0.01] .386 5 .841 1 1.000 -0.01 [-0.02, -0.00] .045 4 .930 0 1.000 0.21 [0.19, 0.22] < .001 2 .986 99 < .001
craving PEV 0.04 [0.03, 0.05] < .001 2 .993 17 < .001 0.04 [0.03, 0.05] < .001 1 .999 15 .008 0.07 [0.05, 0.08] < .001 7 .486 23 < .001
cognitive control ROI 0.02 [-0.00, 0.04] .039 0 1.000 9 .082 0.08 [0.07, 0.10] < .001 0 1.000 12 .007 -0.15 [-0.19, -0.11] < .001 34 < .001 0 1.000
cognitive control PEV -0.07 [-0.07, -0.06] < .001 22 < .001 0 1.000 -0.01 [-0.02, -0.00] .002 2 .991 0 1.000 0.00 [-0.01, 0.02] .371 8 .320 10 .110
craving regulation PEV 0.09 [0.09, 0.11] < .001 0 1.000 17 < .001 0.06 [0.06, 0.08] < .001 0 1.000 10 .031 -0.01 [-0.02, 0.00] .246 1 .993 0 1.000
value ROI 0.05 [0.03, 0.10] < .001 0 1.000 56 < .001 0.04 [0.03, 0.07] < .001 0 1.000 35 < .001 0.07 [0.05, 0.09] < .001 0 1.000 19 < .001
value PEV 0.02 [0.01, 0.04] < .001 7 .564 6 .742 0.01 [0.00, 0.02] .012 1 1.000 1 1.000 0.19 [0.17, 0.21] < .001 0 1.000 81 < .001